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Abstract 


Shifted Legendre orthonormal polynomials (SLOPs) are used to approxi- 
mate the numerical solutions of fractional optimal control problems. To 
do so, first, the operational matrix of the Caputo fractional derivative, 
the SLOPs, and Lagrange multipliers are used to convert such problems 
into algebraic equations. Also, the method is proposed for solving multidi- 
mensional problems, and its convergence is proved. This method is tested 
on some nonlinear examples. The results indicate that the technique can 
efficiently solve multidimensional problems. 
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1 Introduction 


For the first time, fractional calculus was introduced in the 17th century. Li- 
ouville, Griinwald, Letnikov, Riemann, and Caputo substantially contributed 
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to the development of its theoretical foundations [6]. They worked on mass 
and heat transfer problems using the terms semi-derivative and semi-integral. 
The first book on fractional calculus was written by Oldham and Spanier 
[27]. Further details on fractional calculus and some of its applications can 
be found in [11, 12, 21, 22]. 


In recent years, the applications of fractional calculus in engineering and 
sciences, including mathematics, fluid dynamics, and physics, have attracted 
considerable attentions. Fractional calculus is used to extend the usual no- 
tions of derivative and integral to ones with real orders and is based on the 
concepts of fractional derivative in the sense of Caputo and fractional integral 
in the sense of Riemann-—Liouville [22, 27]. 


When we use a term involving fractional-order derivative(s) in differen- 
tial equations of optimal control problems, we obtain fractional optimal con- 
trol problems (FOCPs). Many scientific studies confirm the applications of 
FOCPs in mathematics, mechanics, medicine, and engineering [13, 23]. For 
example, such problems have been used to obtain numerical solutions of 
the fractional models of some diseases, such as the fractional-order tumor- 
immune model, HIV epidemic, and the glucose-insulin system [2, 15, 24]. 


Orthonormal polynomials have been applied in various linear and non- 
linear problems, because they can be used to convert these problems into 
easy-to-solve algebraic equations. They have many useful properties that fa- 
cilitate the solution of mathematical problems and provide a way for solving, 
expanding, and interpreting solutions in some types of differential equations 
[1, 5, 10, 12]. 


In this article, we use the SLOPs as the basis functions of the method 
proposed to solve fractional differential equations. The common approach 
adopted in the past studies was to solve the one-dimensional problem. More- 
over, most of the studies like [5, 4, 10], just obtained the error bound of the 
operational matrix in fractional derivatives. Hence, none of them proved the 
convergence of the method under consideration. 


Therefore, we aim to develop the method for multidimensional problems 
in this paper. Moreover, we prove the convergence of the method. The 
outputs reveal that the method is efficient for multidimensional problems. 


We organized the paper as follows. In Section 2, we present the important 
properties of shifted Legendre polynomials, some preliminary definitions from 
fractional calculus, and the operational matrix of fractional derivatives. In 
Section 3, we explain the method and the necessary conditions for the FOCPs. 
Section 4 discusses the convergence of the proposed technique. In Section 5, 
we compare our results with those of the previous researches for nonlinear and 
multidimensional examples. Finally, in Section 6, we present the conclusion. 
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2 Shifted Legendre orthonormal polynomials 


Definition 1. [5] For a function €(t), the Riemann-Liouville fractional in- 
tegral of order a > 0 is defined by 


wees. eed ae eee, a>0, t>0, 
reo = | rah a= 0, (1) 


where = 
T(a) =} gle * de, 
0 


denotes the gamma function. 


Definition 2. [5] For a function €(t), the Caputo fractional derivative of 
order a is defined by 


1 : ie 
DW = raw f (— 2-01 ee) de, n-l<agn, t>0, 


where n is an integer. 
Some properties of these operators can be written as 


D°*c=0, c is a constant, (3) 


n-1 k 
1° (D* E(t) =) — 1) MOG, (4) 
k=0 : 
oo PG+1) jo 
OTe ia -) 
and 
D% (BE(t) +y7(t)) = BD* E(t) +77 D* v(t), (6) 


where 6, 8, and + are scalar coefficients. 


Definition 3. [3] The Legendre polynomial of degree i, p;(z), is defined on 
the interval [—1, 1] by the recurrence relation 
2i+1 


a . 
piti(z) = Faia | zpi(z) — Fy Pi); a 


W 
B 


(7) 


where 
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po(z) = 1, pi(z) =z. (8) 


We obtain the shifted Legendre polynomials p*(t) on [0,1] if we use the 
change of variable z = 2t — 1: 


. +1 . i, 
po(t)=1, — p(t) = 2t—1. (10) 


These polynomials are orthogonal, in the sense that 


wo.) =f Onoa=fRT Ih an 


As shown in [3], if we introduce the SLOPs p;(t) = /2i + 1p? (t), then 


1 . 3 
(4) a caity: 2% 
[ somme={y 45? (12) 


and 


pi(t) = V2i+1 pee : 1 
k=0 
Assume that ¢ is any element of L?/0,1] and 
pu = span{po(t), Pi(t),---,Pu(t)}- (14) 


Now, for any h € pyz, we can write h ~ yy d; pi(t), where the coefficients 
d; are determined as follows: 


1 
a= f h(t) p;(t) dt, i=0,1,...,M. (15) 
0 


We call ¢, € pu the best approximation of ¢ out of pyy whenever 


for all h € pus: |G — Cpll2 < ||¢ — Alle. (16) 
Since ¢, € pm, there exist coefficients c;,2 = 0,1,..., M, such that 
M 
Colt) S> ci Bi(t). (17) 
i=0 


So, the matrix form of ¢,(t) is 
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C(t) ~ FT A(t), (18) 
where 
co Po(t) 
F= . . Bee 7 (19) 
oT pu (t) 


Theorem 1. For the SLOPs vector Ajs(t), the fractional derivative of order 
a, in the sense of Caputo, is defined as follows: 


D* Au(#) = Dia) Au(2). (20) 


Herein, D(,) denotes the (M+1)x(M-+1) operational matrix of the fractional 
derivative, given by 


0 0 0 0) 
Wig ok 0 0 0 tee 0 
() ~ | Wal(n,0) Wa(n,1) Wa(n,2) --» Wa(n,M) |’ 
W.(M,0) W.(M, 1) WM, 2) fai W..(M, M) 
where 


k j bt ta F : 
= 7 Sy ~ (—1)*tstit! (k + i)! (L +3)! 
=/ RPE CRED DD (—-DMTG— e+ DGG Ga 


i=n 1=0 


and rows 0 to n-1 are zero. 


Proof. See [3]. 


3 The numerical method 


To solve the following problem, we use the operational matrix of fractional 
derivatives, the SLOPs and Lagrange multipliers. 


min J = dh f(t, a(t), u(t)) dt, (22) 
D* x(t) = o(t, x(t), u(d)), n-l<acn,té [t,t], (23) 
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D” x(to) =an, k=0,1,...,n—1. (24) 


Here, $(t, x(t), u(t)) = g(t, x(t)) +0(t) u(t), and S is the feasible solution set. 
Also, u(t) and x(t) denote the control and state variables, respectively, u(t) is 
continuous, x(t) is continuously differentiable, g(t, x(t)), f(t, a(t), u(t)), and 
b(t) are smooth functions, b(t) is invertible, f(t, x(t), u(t)) and d(t, x(t), u(t)) 
are convex functions, S' is a convex set, and f(t, x(t), u(t)) is integrable on 
I = [to,ti]. Moreover, f(t,x(t),u(t)) and g(t,(t)) satisfy the Lipschitz 
property. In fact, 


IIf(t,c1(t), ui(4)) — f(t, 22(4), ua(t)) || < L(|]e1(t) — v2(t)|| + llur(t) on 


and 


lIg(t, x1(t)) — g(t, v2(t)) I] < K (llr (t) — 22(6)II), (26) 


where L and K are positive constants. Approximate x(t) by the SLOPs p;(t) 
as 


Eu(t) =C7 An(t), (27) 
where C7 is an unknown scalar coefficient vector given by 
CT= (co ce cm) : (28) 


We defined p;(t) and Aj,;(t) in (10) and (19), respectively. By (27), we can 
rewrite the dynamic constraint (23) as 


ce Da) Am(t) = g(t,07 Am(t)) + b(t) u(t). (29) 
So, we obtain 
u(t) = nD (CT Dia) A(t) — g(t,C* Am(t))). (30) 


Then, we can rewrite the initial conditions (24) in the form 


CT Dixy Am (to) — te = 9, k=0,1,...,2—1. (31) 


Due to (27), (30) and (31), the performance index J can be approximated 
by 


Ju [C7] = i. f(t,Eu(t), D°Za(t))dt + 3 (CT Dery An(to) — tk) Ak, 
k=0 


to 


(32) 


IJNAO, Vol. 12, No. 3 (Special Issue), 2022, pp 513-532 


Using shifted Legendre orthonormal polynomials for solving fractional ... 519 


where 
1 


F(t, Exc(t), D°Enr(t)) = f(t,C7 Am), nO) (CT Day Au (t) — g(t, C7 Am (t)))), (33) 


and A, denotes the Lagrange multiplier, which should be determined [11]. 
The necessary conditions for the optimality of (22) are subject to the 
dynamic constraints (23) and (24) in the form 


OJu 
Oc; 


OJM 


=0; k=Gi.n wi 6 
OX ’ ree | »n ( ) 


=0,1=0,1,...,M, 


We can use any standard iterative method to solve the aforementioned system 
for c;, i=0,1,...,M, and A,, k =0,1,...,n—1. Asa result, we obtain x(t) 
and u(t) as given in (27) and (30), respectively [3]. 


4 Convergence analysis 


The use of SLOPs operates as a proof of convergence in three steps. In the 
first step, we show that the usage is indeed justifiable. In the second step, 
we show that the functional derivative of a shifted Legendre polynomial is a 
proper approximation for the same derivative. In the third step, we indicate 
the difference between the target function for any optimized solution and the 
value of the target function of the shifted Legendre approximation, tends to 
zero as the number of the shifted Legendre orthonormal basis increases. We 
complete these steps by the hypotheses, Lemmas 1 and 2. To find an upper 
bound for the operational matrix errors in fractional derivatives and to prove 
the convergence, we use the following theorems. 


Theorem 2. Let H be a Hilbert space, and let Y be a finite-dimensional 
subspace of H. Also, assume that {y1, y2,..-, ya} is any basis for Y. Given 
any x in H, let yo denotes the unique best approximation of x out of Y. 
Then, 


GE, Yayenn Yt) 


lla — yolld = Ce (35) 
where 
(z,0) (@,y1) +++ (@,ym) 
CO yc a oe ae (36) 
(um, ) (unr, 1) ie luysastian) 
and 
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(yasy1) o°* (Yi, ym) 
G(yi,Yy2,°°° ,yM) = : : : (37) 


(ym. Yi) ++ (Ym, ym) 


Proof. See [5]. 


We show that the upper bound of operational matrix errors in fractional 
derivatives D‘“) can be obtained as 


e% = D@ Ay (t) — D* Au), (38) 
where D® is an approximation of the operator D( and 


a 
ED,o 
ED 
eH = . ‘ (39) 
a 
ED.M 


As mentioned in [18], for each element of 4, an upper bound for the error 
related to D‘~ can be written as follows: 


k 1 
i k+i)! G(t*~1, po(t), -.-. Bar(t)) \ 2 
leB alle < V2k 4 io a | Copiebiteces . ( Git ») , 
i=1 
0<k<M. (40) 


By Theorem 2 and (40), we conclude that ¢% tends to zero as the number of 
the shifted Legendre orthonormal basis increases [5]. 


Lemma 1. Let x(t) be a continuously differentiable function, and let %y,(t) 
denote the approximation of x(t) by the SLOPs. Then, 


c(t) -Zu(t)|| 30 as M— co. (41) 


Proof. See [15]. 


Lemma 2. For x(t) and %,,(t) as in Lemma 1, when M —> oo, 


|| D® x(t) — D° Ey, (t)|| 0, (42) 
|D* Za¢(to) —2z|=0, k=0,1,...,n—1, (43) 
Z(t) — fm (t)|| > 0. (44) 


Proof. See [5]. 
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We define J1 [C7] as follows: 


JIC? | =i f(t2(t), a (Dia) x(t) — g(t, x(t))) dt 


to b(t) 
+ 5 (Dia) x(to) -— Br) Ak: (45) 
k=0 


Theorem 3. Consider problems (22)—(24), and let x*(t) be an optimal so- 
lution of min J1[C7]. Then, 


| Ju [C7] — J1[C7]| +0 as Mo. (46) 


Proof. Using (27) and (30) we obtain 


ty 
| Jur [C7] — J1 [C7] =|/ f(t,CT Am(®), eh: (CT Dia A(t) — g(t, CT Am(2)))) at 

to b(t) 
n-1 

+ - (c* D xy Su (to) — te) Ak 
k=0 

ty 

= - f(t, x*(t), nD (Dia) x*(t) = g(t, x*(t)))) dt 
n-1 

— So (Bas 2 to) — ee) Ax |: 
k=0 


According to (24), (31), and Lemmas 1 and 2, we know that 


n—-1 
(Cc? Di) Awm(to) — TK)AK =0 
k=0 
and that pee (Do) x* (to) — TK)AK = 0. So, 
| Jur [C7] — J1[C7]| 
ti 
to 
1 
= Pate x(t), o(t) (D(a) x(t) co, g(t, x(t)))))dt| 
We know that f satisfies the Lipschitz condition. Therefore, 
| Jur [C7] — J1[C7]| 
ty 
< f° (E(NCT Asstt) — 200) 
to 


+|| a (C* Dea) Am(t) — g(t,C* A(t) — Dia 2(t) + g(t, 2(4)) ||) at. 
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By the Schwartz inequality and separating integrals, we obtain 
| Jur [C7] — J1[C7]| 


< Lf (Ic? aut ) — a(t))|I) at 


aa fw CP Dea A(t) ~ Dea 2(6)|) a 


mal. (Iig(t,2(4)) — g(t, 07 Aar(t))|)) at 


We write the upper bounds of integrals and note that g satisfies the Lipschitz 
condition. Then, 


| Iu [C7] — J1[C*]| < L(t — to) (IIC* Ans(t) — x(6)I| 


‘ “ba Da) Au(t) — Dea) #(4))|I 


K (ti — to) T 

a I|z(t) -— C° Ane(E)II- 
Dt) 

If M — oo, then Lemma 1 shows that the first and third terms tend to 

zero. Also, the second term tends to zero by Lemma 2. Consequently, 

Ju [CT] > J1[CT]. 


Through Theorem 3, we observed that the difference between the value 
of the target function for any optimized solution of min J1[C7] and that 
of the target function for the approximate value of Legendre tends to zero 
as M — oo. Having (27)-(32) in mind, min J1 [C7] is equivalent to (22). 
Hence, the difference between the value of target function (22) and that of 
the Legendre approximate target function tends to zero. 


5 Numerical experiments 


In this section, we prove the accuracy of the proposed technique by providing 
some examples and then comparing our achievements with the numerical 
results obtained in other papers by the computer with Intel Core i7 CPU up 
to 3.5 GHz, RAM 12GB, and the codes written with Wolfram Mathematica 
11. 


Example 1. Consider the problem 


min T= f (et) — 8) + (u(t) + - 


subject to dynamic constraints 
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D*) «(t) = # x(t) + u(t), (48) 


x(0) = #(0) =0. (49) 
Due to (48), we obtain u(t) and rewrite (47) as 


u(t) = D™ z(t) —#z(2), 


min J -| (ee Am(t) - t)? 


+ (DOT Ay()-@CT Au) +4 —-— 


+ (CT Do) Am(to) — x(0)) Ao + (C* Day Am (to) — &(0))A1. 


£9), 
The functional J is minimized by x*(t) = ¢? and u*(t) = YES) — t*, with 
10 
minimum equal to zero. Table 2 presents the approximate values of J, which 
are obtained by the proposed method and the methods utilized in [21, 3], 
with different values of M. As the results indicate, our approach is better 


than the ones used in [21, 3]. 


Table 1: Approximations of J with different values of 


M | The method | The method used in [21] | The method used in [3] 
4 | 1.66202 x 10-® 6.07530 x 10-6 4.76932 x 10-6 
6 | 2.44576 x 10-7 5.91532 x 10-7 5.37825 x 10-7 
8 | 5.90947 x 10-8 1.21966 x 10-7 1.06099 x 10-7 
9 | 3.26447 x 10-8 7.03371 x 10-8 5.44304 x 10-8 


Table 3 presents the absolute values of errors for the control and state vari- 
ables for various values of t. Also, in Figure 6, the approximate and exact 
values of the control and state variables are plotted for M = 6. 
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Table 2: Absolute errors of x 


Naseri, Heydari and Bagherzadeh 


(t) and u(t) at M =6 


t | I’ @®-2O) 


ju* (t) — ult)| 


0.1 | 1.60241 x 1077 
0.2 | 2.35607 x 1077 
0.3 | 9.96796 x 107-8 
0.4 | 6.68032 x 10-8 
0.5 | 7.86075 x 1078 
0.6 | 9.06389 x 1078 
0.7 | 2.84397 x 107-7 
0.8 | 2.78471 x 1077 
0.9 | 3.55721 x 1078 


1.72334 x 107° 
4.57424 x 1074 
2.85637 x 1074 
2.89849 x 10-4 
1.79588 x 1074 
2.80773 x 1074 
1.15197 x 1074 
2.69036 x 10-4 
2.73064 x 1074 


value 


a 
mie xO 95 


0.2 04 0.6 08 1.0 


0.2 04 06 0.8 


1.0 


u(t) 


Bititg u(t) 


Figure 1: Approximate and exact values of the control and state variables 


for M=6 


Example 2. Consider the two-dimensional problem 


. _ . o 2 ak (ite 3)2 
mind = f ((x1(t) Cy + Get) =f) 


subject to dynamic constraints 


p} x1(t) = 3ur(t) — 3t a1 


M4) a4 T3) jo.9)2 
+ (ur) -# ran’ rag)’ 
T(4) ioe 
+ (u(t) — t+ IF(2.9)' °)?) dt, 


(t) + ? xo(t) — ua(t), 


D"* go(t) = —2 ua(t) + (24 — 1) a(t) + E21 (4), 


and 
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By (51) and (52), we obtain ui(t) and u2(t) as follows: 
“¢ D*1 a (0) —3t?x1(t) + t? x(t) 
i ( ae - lee — 1) a(t) +tar(t) ) 
y(t) = Cy Am(t), CT = (e10 e11 --- crm), 
aro(t) = CF An(t), Cle (c20 C21 “++ Com) 5 


and rewrite (50) as 
min J = [ (( (cP Am(t Pye + (GE Awm(t) -#/ 
1 
+ (5(DY Cf Aw(t) +32 (Cf Au(®) - 2 (Cf Au(®)) 
1 


— =(D™ CF A(t) - (24? - 1) (CP Am()) - t(CT Am(t)) - 4 


M4) ao FE) 
61(2.9) 31°(1.9) 


— (24? — 1) (OF Am(®) -t (CT Am(t))) - 8 + 


on) 


T 


eae ( ou Ce Am(t) 
T(4) 1.9)2 
6r( 9.9) ¢ °)") dt 
+ (Cf Doo) Aw (to) — £1(0))Ao + (cf Dp (1) Awm(to) — %41 (0)) Aq 
+ (CZ De Am(to) — £2(0))Ao + (CZ Dry Am (to) — 22(0)) Ar 


The functions x(t) = ¢?,a%(t) = # and ux(t) = tf — DO 9 4+ 


6T(2.9) 
se 09 us(t) =t — BT 8) t’° minimize the functional J, and the mini- 


mum wale is zero. In Table 4, we present the approximate values of J with 
different values of M. 


Table 3: Approximate values of J with different values of M 


i 

2.39801 x 1077 
3.03043 x 107-8 
6.97336 x 10-9 
6.97321 x 107-9 


cO| oo] or) A] = 


Table 4 presents the absolute values of errors for the state and control vari- 
ables for various values of t. 
Also, in Figures 2 and 3, the approximate and exact values of the state and 
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Table 4: Absolute errors of 71(t), x2 
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t), ui(t), and ua(t) at M =6 


[23 (t) — xo(t) 


Jur (t) — wi (t)| 


[u3(t) — ua(t)| 


t | |zi@)—2()| 
0.1 | 7.19262 x 10-7 
0.2 | 1.0357 x 107° 
0.3 | 3.70976 x 107-7 
0.4 | 4.54804 x 10-7 
0.5 | 5.92208 x 10-7 
0.6 | 9.51419 x 1078 
0.7 | 9.14941 x 107-7 
0.8 | 8.57316 x 107-7 
0.9 | 2.67307 x 107-7 


| 
1.74666 x 1077 
2.48769 x 10-7 
7.82014 x 107-8 
1.37132 x 1077 
1.84041 x 1077 
1.81842 x 1078 
2.02377 x 10-7 
2.32216 x 10-7 
2.16467 x 1079 


6.4603 x 10-5 
1.60228 x 1074 
1.03983 x 1074 
1.05124 x 10-4 
6.48507 x 1075 
1.04023 x 1074 
3.7991 x 1075 
9.89654 x 10-5 
9.10531 x 1075 


9.51622 x 10-6 
2.89678 x 10-° 
1.19302 x 107° 
1.82481 x 1075 
8.03613 x 10~° 
1.65065 x 107° 
6.07151 x 10-6 
1.59221 x 107° 
1.77034 x 1075 


control variables are plotted at M =6. 


value value 

1.0 1.0 

08 08 

06 06 

” X(t) 04 X,(t) 
0.2 minim xf (t) 02 minim x, (0) 


a a «6 2 we” 02 04 06 O08 10 


Figure 2: Approximate and exact values of the state variable at M = 6 


I 

12 value 
. . - Time 
10 06 Os 10 
08 -0.2 
06 u,( 
H -04 * 

04 u,() mite uy(t) 


nitiim u,(t) 96 


Time - 
10 08 


0.2 0.4 0.6 08 


Figure 3: Approximate and exact values of the control variable at M = 6 


We can apply this method to another category of problems. In fact, if in 
problems (22)—(24), we replace (23) by 


yp D® x(t) + Wa(t) = g(t, x(t) +0(t) u(t), (55) 
n-l<acn, b(t) F 0O,te€ [to, t1], 


IJNAO, Vol. 12, No. 3 (Special Issue), 2022, pp 513-532 


Using shifted Legendre orthonormal polynomials for solving fractional ... 527 


then the method still converges according to (44), where y and w are scalar 
coefficients. Let us present one example of this form. 


Example 3. Recall from [28] the problem 


min J = [ (u(t) — x(t))? dt, (56) 
subject to dynamic constraints 
aan 
&(t) + D°x(t) = u(t) — a(t) + Tara + . (57) 
and 
x(0) = 0. (58) 
By (57), we can find u(t): 
7 6 to+2 ’ 
u(t) = &(t) + D°x(t) + a(t) Ta+3) t, 
min J = i (CT Am(t) + D*(CT Am(8)) - ber 3)? dt 
0 T'(a+3) 
+ (CT Do) Am(to) — 2(0))Xo- 
The functions 2*(t) = eS and u*(t) = SS minimize the functional 


J, and the minimum value is zero. In Table 5, we present the approximate 
values of J with different values of M/. 


Table 5: Approximate values of J at a = 0.9 with different values of M 


M J 

4 | 2.32302 x 10-7 
6 | 2.32786 x 10-% 
8 | 2.98816 x 10~!? 


Table 6 presents the absolute values of errors for the control and state vari- 
ables for various values of t. 

Also, in Figure 3, the approximate and exact values of the control and state 
variables are plotted for M = 6. Tables 3 and 8 present the maximum errors 
of u(t) and x(t) with different values of M. 

Also, in Figure 5, the control and state variables are plotted for M = 5 and 
different values of a. 
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Table 6: Absolute errors of x(t) and u(t) at M =6 


|" (t) — x(t)| 


ju* (t) — ult)| 


3.22688 x 10-7 
4.89573 x 1077 
5.31838 x 1077 
6.51328 x 1077 
1.48297 x 1077 
6.3336 x 1077 
1.34478 x 107-7 
5.49314 x 1077 
1.0371 x 107-7 


2.3951 x 1075 
1.18457 x 1075 
1.52362 x 107° 
5.73914 x 10-6 
1.58438 x 1075 
2.83551 x 107° 
1.45402 x 1075 
7.44278 x 10-6 
1.81787 x 1075 


t 
0.1 
0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 
0.9 
value 

0.30 

0.25 

0.20 

0.15 

0.10 


0.2 0.4 06 


x(t) 
sini x*() 


08 


u(t) 


Bittim u(t) 


. , - Time 
02 04 06 O8 0 


Figure 4: Approximate and exact values of the state and control variables at 


M=6 
Table 7: Maximum errors of x(t) and u(t) at M =3. 
M=3 Maximum errors of x(t) | Maximum errors of u(t) 
The method 2.36519 x 10~° 2.30757 x 10-? 
Algorithm 1 in [28] 8.8025 x 10-8 8.8025 x 10-3 
Algorithm 2 in [28] 5.1966 x 1073 4.3260 x 10~? 


Table 8: Maximum errors of x(t) and u(t) at M =5. 


M=5 Maximum errors of x(t) | Maximum errors of u(t) 
Our method 2.21121 x 10-° 4.7773 x 10-4 
Algorithm 1 in [28] 1.0903 x 10-7 1.0903 x 10-7 
Algorithm 2 in [28] A5321 % 10-° 6.3134 x 10-4 
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Figure 5: Control and state variables for M = 5 and different values of a 


6 Conclusion 


In this paper, we applied a numerical method to solve a class of fractional 
optimal control problems. We used the SLOPs and the operational matrix of 
fractional derivatives. Then, we used the Newton iterative technique to solve 
these problems. We obtained the error bound of the operational matrix in 
fractional derivatives and proved the convergence of the method. We focused 
on multidimensional problems, which have never been solved by this tech- 
nique. To show the efficiency of the method for multidimensional problems, 
we provided some nonlinear examples. Comparison of our results with those 
obtained by other techniques in previous studies revealed the accuracy of the 
proposed technique for nonlinear and multidimensional problems. 
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